//!
//!       Staggered Mesh for u-vel and v-vel               Notation
//!      
//!       0   1       2       3       4   5             
//!       |   |       |       |       |   |            o Center of volumes    
//!
//!   5-  #---^--->---^--->---^--->---^---#  -4        > u-velocity 
//!       |       |       |       |       |   
//!   4-  >   o   >   o   >   o   >   o   >            ^ v-velocity 
//!       |       |       |       |       |            
//!       ^---^---+---^---+---^---+---^---^  -3        + nodes
//!       |       |       |       |       |         
//!   3-  >   o   >   o   >   o   >   o   >            # corners, here >, ^
//!       |       |       |       |       |              and o are defined
//!       ^---^---+---^---+---^---+---^---^  -2 
//!       |       |       |       |       |            
//!   2-  >   o   >   o   >   o   >   o   >            
//!       |       |       |       |       |            
//!       ^---^---+---^---+---^---+---^---^  -1 
//!       |       |       |       |       |                 
//!   1-  >   o   >   o   >   o   >   o   >                 
//!       |       |       |       |       |                 
//!   0-  #---^--->---^--->---^--->---^---#  -0 
//!
//!      |       |       |       |       |
//!      0       1       2       3       4                  
//!
//!
//!         Indexation of u (>), v (^) and other variables (o).
//!                                                  
//!                      v(i,j) 
//!                  |     n     |               
//!                --+-----^-----+--             
//!                  |           |               
//!                  |           |               
//!     u(i-1,j) = w >     o     > e = u(i,j)
//!                  |   (i,j)   |
//!                  |           |
//!                --+-----^-----+--
//!                  |     s     |
//!                     v(i,j-1)
//!

//
//---------------------  1D  ---------------------
//

namespace Tuna {

template<class Tprec, int Dim>
inline bool Upwind_CoDiS<Tprec, Dim>::calcCoefficients1D() 
{
    prec_t G_dx = Gamma / dx;
    prec_t dx_dt = dx / dt;
    prec_t ce, cw;
    aE = 0.0; aW = 0.0; aP = 0.0; 
    sp = 0.0;

    for (int i =  bi; i <= ei; ++i)
    {
	ce = u(i  );
	cw = u(i-1);

// The statements:
//	    if ( c? > 0 ) c? = c?; 
//	    else          c? = 0.0;
// Were substituted by:
//	    if ( c? <= 0 ) c? = 0.0; 
// where ? = w and s

	    if ( ce > 0 ) ce = 0.0; 
	    else          ce = -ce;
	    if ( cw <= 0 ) cw = 0.0; 

	    aE (i) = (G_dx + ce);
	    aW (i) = (G_dx + cw);
	    aP (i) = aE (i) + aW (i) + dx_dt;
//		+ (ce - cw) + (cn - cs);
	    sp (i) = phi_0(i) * dx_dt ;
	}    
    applyBoundaryConditions1D();
    return 0;
}

//
//---------------------  2D  ---------------------
//
template<class Tprec, int Dim>
inline bool Upwind_CoDiS<Tprec, Dim>::calcCoefficients2D() 
{
    prec_t Gdy_dx = Gamma * dy / dx;
    prec_t Gdx_dy = Gamma * dx / dy;
    prec_t dxy_dt = dx * dy / dt;
    prec_t ce, cw, cn, cs;
    aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aP = 0.0; 
    sp = 0.0;

    for (int i =  bi; i <= ei; ++i)
	for (int j = bj; j <= ej; ++j)
	{
	    ce = u(i  , j) * dy;
	    cw = u(i-1, j) * dy;
	    cn = v(i, j  ) * dx;
	    cs = v(i, j-1) * dx;

// The statements:
//	    if ( c? > 0 ) c? = c?; 
//	    else          c? = 0.0;
// Were substituted by:
//	    if ( c? <= 0 ) c? = 0.0; 
// where ? = w and s

	    if ( ce > 0 ) ce = 0.0; 
	    else          ce = -ce;
	    if ( cw <= 0 ) cw = 0.0; 

	    if ( cn > 0 ) cn = 0.0;
	    else          cn = -cn;
	    if ( cs <= 0 ) cs = 0.0;
	    
	    aE (i,j) = (Gdy_dx + ce);
	    aW (i,j) = (Gdy_dx + cw);
	    aN (i,j) = (Gdx_dy + cn);
	    aS (i,j) = (Gdx_dy + cs);
	    aP (i,j) = aE (i,j) + aW (i,j) + aN (i,j) + aS (i,j); 
	      //+ dxy_dt;
//		+ (ce - cw) + (cn - cs);
//	    sp (i,j) = phi_0(i,j) * dxy_dt ;
	}    
    applyBoundaryConditions2D();
    return 0;
}

//
//---------------------  3D  ---------------------
//
template<class Tprec, int Dim>
inline bool Upwind_CoDiS<Tprec, Dim>::calcCoefficients3D()
{    
    prec_t dyz = dy * dz, dyz_dx = Gamma * dyz / dx;
    prec_t dxz = dx * dz, dxz_dy = Gamma * dxz / dy;
    prec_t dxy = dx * dy, dxy_dz = Gamma * dxy / dz;
    prec_t dxyz_dt = dx * dy * dz / dt;
    prec_t ce, cw, cn, cs, cf, cb;
    aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aF = 0.0; aB = 0.0; aP = 0.0; 
    sp = 0.0;

    for (int k = bk; k <= ek; ++k)
	for (int i =  bi; i <= ei; ++i)
	    for (int j = bj; j <= ej; ++j)
	    {
		ce = u(i  , j, k) * dyz;
		cw = u(i-1, j, k) * dyz;
		cn = v(i, j  , k) * dxz;
		cs = v(i, j-1, k) * dxz;
		cf = w(i, j, k) * dxy;
		cb = w(i, j, k-1) * dxy;

// The statements:
//	    if ( c? > 0 ) c? = c?; 
//	    else          c? = 0.0;
// Was substituted by:
//	    if ( c? <= 0 ) c? = 0.0; 
// where ? = w, s and b		

		if ( ce > 0 ) ce = 0.0; 
		else          ce = -ce;
		if ( cw <= 0 ) cw = 0.0;

		if ( cn > 0 ) cn = 0.0;
		else          cn = -cn;
		if ( cs <= 0 ) cs = 0.0; 

		if ( cf > 0 ) cf = 0.0;
		else          cf = -cf;
		if ( cb <= 0 ) cb = 0.0; 
		
		aE (i,j,k) = (dyz_dx + ce);
		aW (i,j,k) = (dyz_dx + cw);
		aN (i,j,k) = (dxz_dy + cn);
		aS (i,j,k) = (dxz_dy + cs);
		aF (i,j,k) = (dxy_dz + cf);
		aB (i,j,k) = (dxy_dz + cb);
		aP (i,j,k) = aE (i,j,k) + aW (i,j,k) + aN (i,j,k) + aS (i,j,k) +
		    aF (i,j,k) + aB (i,j,k) + dxyz_dt;
//		+ (ce - cw) + (cn - cs) + (cn - cs);
		sp (i,j,k) = phi_0(i,j,k) * dxyz_dt ;
	    }
    applyBoundaryConditions3D();
    return 0;
}


} // Tuna namespace
